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Relaxation Phenomena in a System of Two Harmonic Oscillators 

Antonia ChimonidoJ*) and E. C. G. Sudarshan 
The University of Texas at Austin, Center for Complex Quantum Systems, 1 University Station C1600, Austin TX 78712 

(Dated: March 7, 2008) 

We study the process by which quantum correlations are created when an interaction Hamiltonian 
is repeatedly applied to a system of two harmonic oscillators for some characteristic time interval. 
We show that, for the case where the oscillator frequencies are equal, the initial Maxwell-Boltzmann 
distributions of the uncoupled parts evolve to a new Maxwell-Boltzmann distribution through a 
series of transient Maxwell-Boltzmann distributions, or quasi-stationary, non-equilibrium states. 
Further, we discuss why the equilibrium reached when the two oscillator frequencies are unequal, is 
not a thermal one. All the calculations are exact and the results are obtained through an iterative 
process, without using perturbation theory. 



I. INTRODUCTION 

n 

The interaction between two isolated systems can be a resource for quantum information processing [lj. However, 
the interaction between an isolated system and the environment surrounding it can lead to decoherence, an undesirable 
loss of information that was initially available. One way or the other, any interaction between two initially uncoupled 
subsystems leads to the exchange of quantities such as purity or polarization 2|, or for thermodynamical systems, 
temperature. To understand and to control this interaction process, as well as to prevent decoherence, we need to 
study the mechanism by which this exchange occurs. There has been a lot of interest and work in this problem 
31 \M> [Si LDj lS| • Most of this work involves the couplin g be tween two (or more) two-level quantum systems, or qubits, 



the fundamental units used in quantum computing 9 
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12|. In this paper, we will instead concentrate on the 



interaction between two harmonic oscillators. More specifically, we are interested in calculating the time evolution 
of each oscillator separately, after an external interaction Hamiltonian has been repeatedly applied to the bipartite 
system for enough successive characteristic time intervals, until a new equilibrium has been reached. 

Consider a bipartite system pi2(0) composed of two initially uncoupled subsystems pi(0) and P2(0). If an interaction 
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Hamiltonian is applied to this system for a characteristic time interval r, the two subsystems will interact and evolve 
to two new states Pi(t) and P2{t). In this paper, we will follow a program in which our system is "refreshed" after 
each characteristic time interval r, and the same interaction Hamiltonian is repeatedly applied to it until the new 
equilibrium is reached. The "refresh" procedure is a very important ingredient in our model and we will spend some 
time explaining it in detail in Sees. II and III. 

This article is structured as follows: in Sec. II, we develop a general program for the relaxation of a bipartite system 
by introducing a periodically repeated interaction scheme. In Sec. Ill, we apply this program to a system of two 
initially uncoupled quantum harmonic oscillators. In Sec. IV, we discuss the physical interpretation and consequences 
which arise from solving this problem and present the conditions needed for a new equilibrium to be attained. We 
explore specific cases in some detail and present our results graphically. Finally, in Sec. V, we discuss the implications 
of our results and suggest possible future directions to be explored. 



II. GENERAL PROGRAM FOR THE EVOLUTION OF A BIPARTITE SYSTEM 



To quantitatively describe physical situations with pure as well as mixed ensembles, we use the density matrix 
formalism, introduced independently by Landau and von Neumann in the 1920's 
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14j . Density matrices not only 



portray the probabilistic nature of a quantum system, but they also contain all the physically significant information 
we can possibly obtain about the ensemble in question. The time evolution of a closed quantum system p(t) governed 
by a time-independent Hamiltonian H is represented by: 

P (t) = U{t,t Q )p{t Q )fr{t,to), 

where U(t, to) is the unitary time evolution operator given by U(t,t ) = er tHt / h . Quantum correlations are generated 
when two such systems are forced to interact through some interaction Hamiltonian Hi nt . 

Consider a bipartite system ^12(0) composed of two initially uncoupled subsystems pi(0) and P2(0). Although the 
evolution of pi2(0) is unitary, the reduced evolution of pi(0) or p2(0) is in general not unitary [15j . Knowing the form 
of p\ at any time t\ is not sufficient to predict its form at a future time ti. So to extract the reduced evolution of 
subsystem 1, we need to first calculate the time evolution of pn{$) and eliminate the effect of the unwanted part by 
taking the partial trace of the evolved complete system over subsystem 2. Namely: 

Pi(t)=Tr 2 [U(t)p 12 (0)U f (t)l 

where pi2(0) = pi(0) <8> ^2(0) and U(t) — eT %m l n , with H being the total Hamiltonian of the composed system, 



including H 



int • 



16] , in which a constant interaction Hamiltonian 



In our work, we will follow a scheme similar to one proposed by Rau 
Hi n t is applied to a bipartite system pi2(0) = Pi(0) ® P2(0) composed of two initially uncoupled harmonic oscillators 
with density matrices pi{0) and p2(0), and with p2 kept in an inexhaustible temperature bath. In Rau's model, 
emphasis is given to studying how system 1 relaxes when it is in contact with system 2. The interaction Hamiltonian 
is applied to the composite system pi2(0) for a constant characteristic time interval t, interrupted, and then reapplied 
for another time interval r to the time-evolved state Pi2(t) = Pi(t) <g) P2,{t), where pi(r) = Tr2[^7(r)/9i2(0)C/^(r)], and 
P2{t) = /02(O). As mentioned above, system 2 is assumed to be in contact with an inexhaustible temperature bath, 
or reservoir, so that any change it may undergo during the time interval r is negligible. In other words, at t > 0, 
a general, constant interaction Hamiltonian is applied for a characteristic time interval r, the system is allowed to 
interact, it is then "refreshed" , and the process is repeated. The "interaction" process is defined as follows: the two 
quantum systems pi(0) and p2(0) interact with one another through the interaction Hamiltonian Hi nt during a time 
interval < t < r, where r is fixed. In Rau's model, the "refresh" process is the extraction of the reduced time-evolved 
quantum part p\ (t) from the complete system p\2 (t) and the resetting of p2 (j) back to its initial form. Finally, the 
"repeat" process is the application of the same interaction Hamiltonian to the new state p\2 (r) = p\ (r) ® p2 (t) (with 
P2{t) = /02(O)) for another time interval r < t < It. There are several areas in physics where such a model is a good 
approximation. For example, in spintronics, in the case of paramagnetic relaxation, a system of spins interacts with a 



lattice, and the information is encoded in the spin states of the conduction electrons 17[ • In this case, the conduction 
electron takes the role of system 1, and the atoms in the lattice take the role of system 2, or the environment. As 
the conduction electron passes through the lattice, it interacts with a lattice atom, and as a consequence, its state is 
changed. This new state will now interact with another atom in the lattice, identical to the first one. 

The following assumptions are made in the formulation of generating relaxation as described in Rau's model: First, 
it is assumed that there is statistical independence between the two density matrices of subsystems 1 and 2 at t = 0, 
or in other words, the density matrix of the combined system can be factorized into density matrices of the component 
subsystems at t = 0, i.e., puiO) — Pi(0) <8> /?2(0). Second, the relaxation time of the system of interest, subsystem 1, 
is much larger than the characteristic interval r. In other words, the system is reset before the time interval specified 
by the inverse of the interaction energy. Third, time averages are never taken, but after every characteristic time 
interval r, partial trace with respect to the reservoir is taken. Finally, system 2 is assumed to be in an inexhaustible 
temperature bath so that any change in it during the interaction process is negligible. At each time interval r, the 



system goes back to its original state at t = while system 1, the system of interest, changes its state. The mechanism 
of refreshing introduced in this model annuls the correlation between the two systems 1 and 2 at every characteristic 
time interval r and replaces the time-evolved state of subsystem 2 with an identical copy of its original state. 

Implicit in this model is the assumption, invoked by Pauli |18j, that the occupation numbers, which correspond 
to the diagonal elements of the density matrix, remain good quantum numbers at all times. In other words, it is 
assumed that, during the course of time, off-diagonal terms become unimportant. Pauli eliminated them by invoking 
a "random- phase approximation" at all times. This approximation is a serious ingredient in the analysis of this 
problem and we believe it is worthwhile to analyze its use at this point. To begin, we stress that the disappearance 
of interferences between a system and its environment is an irreversible process. Irreversible, or energy-dissipating 
processes always involve transitions between quantum states. Such processes are described, at the simplest level, by 
master or rate equations. The Pauli master equation [18j . is the most commonly used model of irreversible processes 
in simple quantum systems. It can be derived from elementary quantum mechanics and by neglecting off-diagonal 
terms. Among other (inessential to our discussion) approximations, the derivation of the Pauli master equation suffers 
from the restriction that one has to discard any built-up of phase relations by invoking a repeated "random-phase 
approximation" at a series of times at microscopically small intervals. This restriction may find itsjustification in Van 
Kampen's analysis of the problem in deriving the master equation from quantum mechanics 20| . Specifically, Van 



Kampen observed that, by writing a master equation, we only intend to get statements on macroscopically observable 
properties of (statistically) large systems. In trying to derive a master equation from quantum mechanics, we must 
therefore first construct, following Pauli, a suitable coarse-graining of phase-space in such a way that the quantities 
of the statistical theory are only those that can be measured macroscopically. Van Kampen's derivation of the master 
equation highlights quite explicitly the inherent difficulties of non-equilibrium statistical mechanics. To derive a 
master equation, we must postulate a number of (justifiable) mathematical assumptions, unfortunately in most cases 
without being able to give the explicit criteria on the microscopic dynamics of the system for these assumptions to 
be valid. The general attitude is that because many large systems evolve smoothly on a macroscopic time scale, 
microscopic details are most likely not important, and must therefore be suppressed. Further, it is assumed that the 
evolution of the reduced density matrix is a Markoff process. This translates to invoking a repeated random- phase 
approximation, i.e. neglecting or suppressing any dynamical built-up of phases as time evolves. For example, in the 
case of photon scattering, interference terms connecting different positions become unobservable at the macroscopic 
body itself, though still existing in the whole system. This was first discussed by von Neumann in his theory of the 



measurement process 



211 ] . Another possible physical mechanism which may aid in explaining the validity of Pauli's 



22], |23|, that the 



approximation is the interaction of the system with its natural environment. It has been shown 
interaction between the system and its environment causes some interference terms to become unobservable. 

In this paper, we will follow a scheme that deviates slightly from the one proposed by Rau. In our model, we 
again start with two initially uncoupled subsystems described by density matrices pi(0) and P2(0). The crucial 
difference between the program analyzed by Rau and what we propose here, is that system 2 is no longer kept in 
an inexhaustible temperature bath. After a constant interaction Hamiltonian is applied to the bipartite system for 
a characteristic time interval r, the system is "refreshed" , but the "refresh" procedure is no longer the same as that 
described above. Instead of concentrating our attention to studying how system 1 relaxes when in contact with system 
2 like in Rau's model, we will now look at how each system relaxes when in contact with the other, once H int is 
applied to the combined system. System 2 will no longer act as the reservoir which goes back to its initial state after 
every r. Instead, both subsystems will be equally important and both will undergo a change in their states during 
the interaction process. We will no longer be concentrating on a system in contact with an unchanging environment, 
but in two systems of interest in contact with each other. 

The "refresh" procedure as defined in our model can be written as a transformation p\2 — > Tri\pv^\ ®Tr\[pi2\. As 
an aside that will become transparent later on, we state at this point, that even though the "refresh" procedure is 
expressed as the tensor product between the time-evolved reduced subsystem states obtained by performing partial 
traces over each unwanted part separately, the operation of taking the partial trace is performed only once, and over 
only one of the oscillator subsystems. The partial trace operation assumes that we consider the relative phase between 
systems 1 and 2 to be completely randomized and thus we can integrate over it. In addition to this, we assume that 
we know nothing about the occupation numbers of each mode of the second system and hence we can sum over them 
when computing the partial density matrix of system 1 (and vice- versa for system 2). Further, we assume statistical 
independence between the two density matrices p\ and pi at each time r. For the time being, we simply state that 
the "refresh" procedure is composed of nothing but a single spectrum measurement in the number basis of one of the 
oscillators. This single measurement has a twofold effect: it diagonalizes the output global state in both oscillator 
number bases {|ni)} and {|n2)}, and at the same time decouples the two oscillators by allowing us to write the global 
output state as a tensor product of the two reduced, post-measurement states. The physical implementation of our 
model along with this last assumption will become clear in Sec. Ill, where we explicitly calculate the reduced density 
matrix pi(r), and extensively discuss the details of the "refresh" procedure. 



Here, we must emphasize that the "refresh" procedure defined in our model is different from the Pauli "random-phase 
approximation" employed in Rau's model. Formally, the "refresh" procedure is a transformation p i2 — > Tr 2 [pi 2 ] <8> 
Tri[pi2] and appears to be nonlinear with respect to p i2 (this can be easily observed if pi 2 is written as a convex 
sum of pure density matrices). In contrast, the Pauli "random-phase approximation" is a transformation p i2 — > 
S ni n 2 \ n i n 2) {nin2\pi2\n\n2) {nxn 2 \ (with |nin 2 ) being the orthonormal basis of the total free Hamiltonian), that is 
linear with respect to p 12 , and corresponds to vanishing of off-diagonal terms in p 12 . Whereas the partial trace assumes 
that we have no knowledge of the occupation numbers of each mode in the second subsystem and hence we end up 
averaging over those occupation numbers, the "random-phase approximation" considers all phase relationships to be 
random. In contrast to Rau's model, we do not have a Markovian master equation for either of the two oscillators. 
Since, in our case, the second system is not reset back to its initial state after each time interval r, the other system 
keeps a memory of what happened in the previous step. The nonlinearity of the "refresh" procedure may raise 
philosophical questions on the physical realizations associated with it. At this point, it is perhaps appropriate to 
clear any confusions that may arise due to the apparent nonlinearity of the transformation: the "refresh" procedure is 
expressed as the tensor product between the two reduced post-measurement states, written as states which are obtained 
by taking partial traces over the corresponding unwanted parts. However, like mentioned above and discussed further 
in section III of this article, the transformation can be obtained by making a single, spectrum measurement on the 
number basis of just one of the two oscillators. This measurement will not only determine the probability distribution 
of the occupation numbers in oscillator 1, but in addition, that of oscillator 2. We strongly emphasize that the 
operation of taking the partial trace is performed only once. This procedure is no different than that of measuring 
the spin of an entangled electron pair: a measurement performed on one of the electrons not only determines that 
electron's spin, but also the spin of the second electron. Therefore, even though in a sense the "refresh" procedure 
appears to be nonlinear in pi 2 , employing it does not break any physical rules of quantum mechanics. 

That being clarified, the "interact-refresh-repeat" process of our model is described in more detail below: 

Pi2(0) = pi(0) <g> p 2 (0) ^ pi 2 (r) = pi(r) <g> p 2 (r) ^> . . . ^> pi 2 (nr) = pi(nr) <g> p 2 (nr). 
More specifically, 

Pi(0) ® flj(O) ^ Tr 2 [f7p 12 (0)f/t] ® Th[#pi 2 (0)#+] ^ . . . ^ Tr 2 [u Pl2 [(n - 1)t]&] ® Th [u Pl2 [(n - l) T ]fr 
with the time evolution operator given by: 

U( T ) = e- lkT ' h . 



A flow chart of our model follows in Fig. [TJ 



7 




Tr 2 p 12 (x) 



Tr, p 12 (x) 



Repeated n times 



FIG. 1: Flow Chart describing the first step in the iterative process of "interacting", "refreshing", and "repeating". Two initially 
uncoupled subsystems are forced to interact through a constant interaction Hamiltonian Hi n t, the time-evolved combined density 
matrix is calculated, and the system is "refreshed". Statistical independence is assumed at each characteristic time interval r. 
The interaction Hamiltonian is again applied to the new uncoupled states of the two subsystems. The process is periodically 
repeated until a new equilibrium is reached. 

This iteration process is carried out until a new equilibrium is reached. Our aim is to study the large time behavior 
of this periodically repeated interaction scheme and understand its effect on the reduced quantum subsystems pi and 

P2- 

The problem of describing the approach to equilibrium for systems composed of a large number of interacting 
particles is a very fundamental one, and one that has been studied extensively through the years. Ever since the 
classic work of Boltzmann, it has been recognized that not only statistical, but also dynamical considerations play a 
role in the time evolution of these systems. Stochastic processes like the one described above are encountered in many 
contexts in both physics, chemistry, and, more recently, quantum communication. Relevant experimental applications 
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involve the investigation of resonance Raman scattering 2J, |25[ , as well as collisional decoherence during writing and 
reading quantum states [26 1 . 

III. EVOLVING HARMONIC OSCILLATOR SYSTEM 

We are now ready to apply the above recipe to the case of two coupled harmonic oscillators. The system considered 
here is a set of two initially uncoupled harmonic oscillators 1 and 2, described by density matrices pi(0) and P2(0) 
respectively. The two subsystems are initially at equilibrium at temperatures Ti(0) and T" 2 (0). The free and interaction 
Hamiltonians of our system are given by: 

Hi = hujia\a\ 
H-2 = hD 2 a\a2 
Hi n t = hjj\(a[a2 + a\a{), 

where wi t 2 is the oscillation frequency of oscillator 1, 2, and ui is the frequency of the applied interaction Hamiltonian. 
A is the strength of coupling between the two oscillators. Here, a\ and a\ are the raising and lowering operators for 
oscillator 1, and similarly, a\ and a 2 are the raising and lowering operators for oscillator 2. For simplicity, we will set 
H = 1 for the remainder of this paper. 

The initial uncoupled and combined density matrices have the Maxwell-Boltzmann distribution and are given by: 



Pi(0) = 



ZxMo)] 



-wid}di0i(O) 



0o ((\) - p -w2a|a 2 e 2 (o) 

P2[) Z2MO)} 

Pl2lUj z 1 [0 1 (o)]z*Mo)] 

= -[wiaJai0i(O)+W2a^a 2 e2(O)] /i\ 

where we have used the property [djai,ci2^2] = to obtain the last equation. Here, 812 — l/(fcT 12 ), k is the 
Boltzmann constant and Ti j2 is the initial equilibrium temperature of oscillator 1, 2. The quantity Z is a normalization 
constant, the partition function, and is given by Z l = e -E„/(kT) = ^^-^/(fer)^ = r^^-u&l^y Thc timc 
evolution of the complete system is given by: 

= e -i{Hi+H2+H tn t)r pi2 {Q) e +i(Hi+H 2 +H mt )r 



1 



Zi[0i(O)]Z 2 [0 2 (O)] 



x e 



■[wiataiei(0)+a;2ata2e2(0)] 



X g+ ir [ w i a i a i+ w 2d2a2+wA(a{a2+d^ai)] ^) 

We will eventually extract the reduced evolution of one part by taking the partial trace with respect to the other 
part, but first we need to calculate pi 2 {r) explicitly. To proceed, we note that there is a very interesting connection 
between the algebra of angular momentum and the algebra of two uncoupled harmonic oscillators. This connection 
is the well-known Schwinger's Oscillator Model of Angular Momentum [27J. In this model, we define the following 
operators: 



J+ 


— a\a 2 




J_ 


= a\a\ 




Jl 


= 2*- a i a2 + 




k 


1 rt ~ 


- a\a\) 


k 


l r t- 


a 2 a 2 ) 


j 


l 

= 2( a i fll + 


a\a 2 ) ■ 


m 




a\a 2 ) ■ 



2 V 

-(tlx - n 2 ) = J 3 . 

It can be shown that the above operators satisfy angular momentum commutation relations. With the above defini- 
tions, equation © becomes: 

_ ( T \ = \ -i[j(aJi+w 2 ) + J3(wi-W2)+2a)AJi]T 

Pl2[ ' z 1 [e 1 (o)]z 2 [e 2 (o)} 

x e -j[wie 1 (o)+w 2 82(o)]+J3[wie 1 (o)-w 2 e2(o)] 

x e +*[j( w l+"2) + .-?3(wi-W2)+2wAJi]r 

Since j commutes with J\ and J3, we can pull it through and factor it out, obtaining a simplified expression for pi 2 (r): 

mn ( T \ = I p -j[t»i6i{p)+w 3 e2(o)] 

Pl2[ ' Zi[0i(o)]z a [0 a (o)] 

X g — "'[•Mk'i— a) 2 )+2a>AJi] 

x e -J 3 [wiei(o)-w a e 2 (o)] 

X e + iT lh(<^i-^2)+2u)\J 1 ] 
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By making the following definitions: 



a = t(cji — uj 2 ) 

b = 2uj\t 

c = i[w 2 e 2 {0) - u 1 1 {0)], (4) 
we can rewrite equation ^ into its final form as: 

n ,oM - - p -j["i0i(O)+a; 2 2 (O)] -i(aJ 3 +bJi) -icJ 3 +i(a.J 3 +bj 1 ) (k\ 

The reduced time evolution of oscillator 1 after a time r is obtained by taking the partial trace of equation (JSJ over 
oscillator 2. The matrix elements of such an expression are given by: 



A\Pi{ T )\ n i) = ^2Win 2 \pi2(r)\nin 2 ) = S n > ini ^2(n' 1 n 2 \pi 2 (r)\n 1 n 2 ) 

"2 

-j[wi0i(O)+a»2ea(O)] 



no 

1 



Zi[0i(O)]Z 2 [0 2 (O)] 



2j(nin 2 |' 



11 2 



x e ~ l{aj3+bJl) x e~ lcj3 x e +l(Q ' /3+6Jl) |ni«2), (6) 

As promised in section II, we now discuss the physical realizations of the "refresh" procedure. To begin, we note 
that the reduced density matrix pi(r) is diagonal in the number basis {|«i)} of the first harmonic oscillator. This is 
due to the fact that j = N = l/2(o]ai + a 2 a 2 ) is a conserved quantity in our model ([if, N] — 0), and also because 
the initial density matrix pi 2 (0) is diagonal in the number basis {|win 2 )} of the harmonic oscillators. The evolution 
generated by the interaction Hamiltonian cannot change the value of the conserved quantity. Since the only non-zero 
elements in pi 2 (0) are those satisfying ni+n 2 = n^+n^, the only non-zero matrix elements of the time-evolved density 
matrix pi 2 (r) are (nin 2 \pi 2 (T)\n' 1 n 2 ) with n\ +n 2 = n' 1 +n' 2 . In performing the partial trace over the second harmonic 
oscillator, we set n 2 = n' 2 , which means that all the terms of pi 2 (r) contributing to the partial trace must also have 
rii — n[. Like mentioned in Sec. II, the operation of taking the partial trace over the second oscillator is equivalent 
to making a non-selective, or spectrum measurement in the number basis of the first oscillator. It so happens that 
because of the conservation of the global number of particles, this one spectrum measurement on subsystem 1, not 
only determines the probability distribution of the occupations numbers of the first harmonic oscillator, but it also 
automatically determines that of the second oscillator. In addition, it simultaneously diagonalizes the output state in 
both number bases {|^i)} and {|n 2 )} of the two oscillators. Since the output state is now diagonal in both number 
bases {|^i)} and {|n 2 )}, it can be written as the tensor product of the two reduced post-measurement subsystem 
states and statistical independence between the two subsystems can be assumed. Even though the "refresh" procedure 
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is expressed as the tensor product of the two post-measurement states, obtained by taking partial traces over each 
unwanted part, the actual operation of taking the partial trace is performed only once, and over only one of the 
subsystems. A single measurement determines the outcome of not just one, but both oscillator systems. 

We now need to figure out how the exponential operators in equation © act on our two-system state \n\n2)- This 
calculation is greatly simplified if we make a transformation to the Euler-angle form e r l0lJ ^e~ % ^ j2 e -ryj3 , since the 
action of the operators J3 and J 2 on our state is a well-known result. We would like to find expressions for the 
coefficients a, /3, and 7 that appear in the Euler angle form, in terms of the coefficients a, b, and c given in equation 
(HI). To do this, we make use of the well-known identity: 



e m ' a 2 = I cos — — in- a sin — 



where, for our case, 



<t = o\i + a-2] + ask = 2 J, 



(the 171^,3 are the Pauli spin matrices) and 



such that 



b - a 

(f)2 + ( | )2 V- . 



(5 



After some extensive algebraic calculations, the results are: 



:d\ ( b 

a = arctan 1^1+ arctan I — — 



(3 = 2arccos 



.4 



cos (arctan -j) 



' D\ ( B 

7 = arctan [ — ] — arc; an 



C ' 



where 



.4 
B 

C 

D 



cos 
ab 



2d 2 
b 



(I) 

iin 2 {d) sin (|) 



cos(d) sin(d) sin ^ ? 



cos 2 (c!) 



d 2 



■sm 2 (d) 



and 



d = 



2 



With the above transformations, equation ([6]) becomes: 



KlPiWK) = 



- / n 1 ni 



Z 1 [0 1 (O)]Z 2 [0 2 (O)] 

1 



£«-* 



(n 1 +n 2 )[wi9i(0)+w 2 e2(0)] e -^(ni-n2)(a+7) 



x(nm 2 |e I/3j2 |nm 2 ), 



where the J 2 matrix element in the last equation is given by Wigner's Formula in the \j,m) basis 28 1 
0\m'|e-^|j,m> = d« m (/3) 

- £(-d 



\/ (i + m ) ! (i - + - m ') ! 



A" 



(j - ml - K)\K\{j + m - K)\(K + m! - m)\ 



P \ j-\-Tn-\-j—Tn —IK / p\ 2K-\-m—rn / 

x ^cos-J (^sin-j 
Substituting this in equation j7]) and summing over n 2 and K, we obtain: 



KM r )l™i) 



8 n 



1 — COS ( ^ I r - 



Zi[0i(O)]Z 2 [0 2 (O)] 

x>s \ / 1 - sec (§) ei( Q 



£ ) P 4(a+7)-e 



+ 7 )-e 



1 — cos I 5 I < - 



£ ] P i(«+7)-e 



where 



= - Mi(O)+wa0a(O)]. 
Defining a = i£ and 7 = i£ in equation ([5]) and letting <5 = |(£ + £), we finally get: 



(n' x \pi(T)\nx) 



Equation @ can be written in the form: 



Zi[0i(O)]Z 2 [0 2 (O)] 

1 

1 -cos(|)e-(e+*)' 



cos(f )e s+e - 1 



*e-s[ e e+s _ C os(f 



Z[Vi[t)\ 



-uiiniBx(r) 



with 
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and 



z[6 1 {t)] = z l [e 1 {Q))z 2 \e 2 {Q)] 



1 — cos( — )e 



We could have obtained the matrix elements of the reduced density matrix of oscillator 2, by taking the partial trace 
of the time-evolved density matrix of the complete system over oscillator 1, by following a similar procedure as the 
one presented one. 

Equation (|10p is a quasi-stationary, non-equilibrium state. It states that, after the interaction Hamiltonian has been 
applied for a time interval r, the original Maxwell-Boltzmann distribution describing oscillator 1 with temperature 
Ti(0) = 1 /[/c6*i (0)] at t = 0, is replaced by another Maxwell-Boltzmann distribution and a different temperature 
T\(t) = l/[fc^i(r)]. More specifically, the interaction has forced the oscillator out of equilibrium at temperature Ti(0) 
and brought it to a new equilibrium at temperature Ti(r). A similar statement is true of oscillator 2. 

Now that we have the matrix elements of the reduced time-evolved density matrix of oscillator 1 at time r, we 
can repeat our calculations to obtain (?i' 1 |pi(2T)|ni), which will lead to another Maxwell-Boltzmann distribution 
(and another quasi-stationary, non-equilibrium state), with a different temperature 7i(2r) = 1/[/c#i(2t)]. The above 
process can be repeated to obtain (n' 1 |pi(3r)|ni), (n' 1 |pi(4r)|ni), (n' 1 |pi(nr)|ni) with an increasing level of 
difficulty. 



IV. PHYSICAL INTERPRETATION AND CONSEQUENCES 

The question we want to answer now is: When is equilibrium achieved! To answer this question, let's consider the 
system composed of some energy state |ni) in oscillator 1 and some energy state |n2) in oscillator 2. Equilibrium will 
be reached when the net rate of "absorption" equals the net rate of "emission" in and out of our system. "Emission" 
results from transitions such as \rii) — ► \n\ + 1} or — > — 1} in oscillator 1, and \n 2 ) — > \n 2 — 1} or \n 2 ) — * \n 2 + 1} 
in oscillator 2. Likewise, "absorption" processes are caused from transitions such as \ni — 1} — > \n\) or \n\ + 1} — > \ni) 
and |n 2 + 1} — > \n 2 ) or \n 2 — 1) — » \n 2 ) for \n 2 ). 

Recalling that Hint = uiA(a\a 2 + a\a{), we can calculate the probability amplitudes for "emission" out of this 
bipartite system: 

(ni + 1, n 2 - l\u)\(a\a 2 + a\a\)\ni,n 2 ) = u\^n 2 {ni + 1), 

and 

(ni — 1,712 + l\uj\(a\a 2 + a2,di)|7ii, n 2 ) = ni(n 2 + 1). 
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Similarly, the probability amplitudes for "absorption" are given by: 

(ni,n 2 |wA(a{a 2 + a\ai)\m — l> n 2 + l) = u)Xy/m(n2 + 1), 

and 

{n 1 ,n 2 \u>X(a\a 2 + a|oi)|ni + l,n 2 - 1} = ui\y / n 2 (n 1 + 1). 

At equilibrium, we expect that: 

P(rci,n 2 )[n 2 (rci + P{n ll n 2 )[n 1 {n 2 + 1)] = P{n x -l,n 2 + l)[nx{n 2 + 1)] +P(m + l,n 2 - l)[n 2 (m + 1)], (11) 

where P{n\, n 2 ) is the probability to have a particle in energy level \ni) in oscillator 1 and energy level \n 2 ) in oscillator 
2, etc. 

Simplifying equation (|11[) gives: 

P(ni,n 2 )[2nin 2 + rii + n 2 ] = P(«i — 1, n 2 + l)[ni(n 2 + 1)] 

+ P(ni + l,n 2 -l)[n2(ni + l)]. (12) 

We can now explore the consequences of equation (fT2"| . First, we note that the interaction Hamiltonian imposes 
the condition that the sum of the occupation numbers n\ and n 2 in states |ni) and \n 2 ) respectively, is a constant. 
Let's call this constant N, such that n\ + n 2 = N. The case N = is trivial, it represents the vacuum state. For 
N = 1, there are two possible transition states, 1 1 , 0) and |0, 1). Solving equation (jT3J) for these two different cases 
leads to the condition P(1,0) = P(0, 1). This result is not surprising. It simply states that, at equilibrium, the 
probability of a transition to the state 1 1 , 0) is the same as that to the state |0, 1). Similarly, for N — 2, there are 
three possible transition states: |0,2), 1,1), and |2,0). This time, equation (JT3J) yields a more interesting result: 
P(0, 2) = P(l, 1) = P(2, 0). Since the sum of all probabilities must be unity, each of these must be 1/3. Repeating 
this procedure for increasing N, we soon recognize that for any N, there are N + 1 different possible states, and the 
probability of transition to any of them is the same for all and equal to 1/(N +1). 

Making use of the above result and employing equations fTJ) at t — and (fTUI) at equilibrium, we conclude that the 
new density matrix for the combined system will be of the form: 

(ni.nahaOxOlm.na) = J_ e -»>»»M«>) e H*nrf J (oo) i (13) 

Recalling that N = n\ + n 2 and defining v — n\ — n 2 , we can rewrite equation P^|) as: 

,N + v N-v. , . ,N + v N-v. 

( 7T~ ' o ^12(oo) I , ) 
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J_ p -"i(^)M°°)p-^(^)0 2 (oo) 



N + 1 
1 

N+l 



e -^-[uiOi(oo)+ui202(oo)] e — Si (oo) -W202(ao)] 



(14) 



Since, for any AT, the probability of transition at equilibrium to any of the N + l possible states is equal to 1/ (N + 1) , 
equation (|14[) should be independent of v. This is true when [wi Si (oo)—ct> 2 02(00)] _^ ^ which leads to the consequence 
that, at equilibrium, Wi0i(oo) = ^2^2(00), or o;i/[fcTi(oo)] = ^/[fclMoo)]. Note that this result is irrespective of the 
initial temperature of the oscillators. 

In the plots below, the "temperature" (kto/H)T and the quantity kT/huj are plotted as functions of the refreshing 
time intervals tit /to, for the relaxation of the initial Maxwell-Boltzmann distributions of the two oscillators. Here, 
we introduce the time constant to for the simplicity of working in a dimensional system of units. In this system, we 
set k = h = to = 1. The interaction Hamiltonian is applied for a constant time interval r, interrupted, the system is 
refreshed, and the procedure is repeated, until a new equilibrium is reached. In Figs. [2]and[3l the frequencies of both 
harmonic oscillators and that of the interaction Hamiltonian are equal, in Fig. 01 the frequencies of the oscillators are 
equal, but that of the interaction Hamiltonian is not, and, finally, in Figs. [5]and[6l the frequencies of the oscillators 
and the interaction Hamiltonian are all different. The dashed line represents oscillator 1 and the solid one represents 
oscillator 2. 

T 1 





5 10 15 20 25 30 5 1 15 2 2 5 30 

(a) (b) 

FIG. 2: (kto/h)T and kT/hw as functions of the refreshing time intervals nr/to for the relaxation of the two initial Maxwell- 
Boltzmann distributions. The dashed line represents oscillator 1 and the solid one represents oscillator 2. Here, k = K = to = 1, 
(kto/h)Ti(0) — 1, {kto / h)T2(0) — 9, uiito = 1^2*0 = uto = 1, and r/to = 2.7. The frequencies of the oscillators are the same 
and the equilibrium reached is a thermal one. 
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4 6 2 4 

(a) (b) 

FIG. 3: (kto/h)T and kT/hw as functions of the refreshing time intervals nr/to for the relaxation of the two initial Maxwell- 
Boltzmann distributions. The dashed line represents oscillator 1 and the solid one represents oscillator 2. Here, k = h = to = 1, 
(kto/h)Ti(0) — 1, (kto/h)T2(0) = 9, Wito = W2<o = wto — 1; an d f/io = 0.3. Keeping the oscillation frequencies the same as in 
Fig. [2] the characteristic time interval is decreased and the equilibrium is reached much faster. 

T I 





FIG. 4: (kto/h)T and kT/fko as functions of the refreshing time intervals nr/to for the relaxation of the two initial Maxwell- 
Boltzmann distributions. The dashed line represents oscillator 1 and the solid one represents oscillator 2. Here, k = h = to = 1, 
(kto/h)Ti(0) — 1, (kto/h)T2(0) — 9, Wiio = ^2*0 = 1, <^io = 5, and r/to — 2.7. Keeping the characteristic time interval the 
same as in Fig. [2j the interaction frequency is increased and the equilibrium is reached much faster. 

We note some interesting results. The rate with which the equilibrium is reached depends on an interplay between 
the characteristic time interval r and the interaction frequency lo. Comparing Figs. [2] and [5J we observe that, for 
a constant interaction frequency and a shorter characteristic time interval, the harmonic oscillator system reaches 
equilibrium much faster. This can be understood by noting that as r decreases, the period of the applied interaction 
Hamiltonian is decreased, and the oscillators do not have as much time to interact before the interaction is interrupted. 
Decreasing the length of the interaction prevents the oscillators from returning back the temperature they have 
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i 




1.5 2 2.5 3 1 2 3 4 

(a) (b) 

FIG. 5: (kto/h)T and kT/hw as functions of the refreshing time intervals nr/to for the relaxation of the two initial Maxwell- 
Boltzmann distributions. The dashed line represents oscillator 1 and the solid one represents oscillator 2. Here, k = h = to = 1, 
(kto/h)Ti(0) — 2, (kto/KjTzip) — 6, uiito = 1, CJzto = 3, uito = 5, and r/to = 1-5. The oscillator frequencies are unequal and 
the equilibrium reached is not thermal. Figure (a) shows that no single temperature can be assigned to the new equilibrium 
distribution, and Fig.(b) demonstrates that the condition uj\/[kT\(oo)] = ui2/[kT2(oo)] is satisfied. 



IftA'ir.f^v^-.v— — - 
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1 2 3 4 5 6" 1 2 3 4 5 6 

(a) (b) 

FIG. 6: (kto/h)T and kT/hu as functions of the refreshing time intervals nr/to for the relaxation of the two initial Maxwell- 
Boltzmann distributions. The dashed line represents oscillator 1 and the solid one represents oscillator 2. Here, k = h = to = 1, 
(kto/h)Ti(0) — 8, (kto / h)T2(0) = 2, uiito = 10, u>2to = 4, u)to = 5, and r/to = 1.5. Again, the new equilibrium distribution 
does not satisfy Maxwell-Boltzmann statistics and no temperature can be assigned to it. 

exchanged so the equilibrium is achieved faster. The same effect is observed when we compare Figs. [2] and 21 
for a constant characteristic time interval and a higher interaction frequency, the oscillators are forced to exchange 
temperature much faster. Again, the equilibrium is attained more quickly. 

Further, we note that our prediction is verified: Wi#i(oo) = w 2 ^2(°o) ( or - w i/[^7i(oo)] = ^/[fcX^oo)]), at equi- 
librium. We observe that when u/i = u>2, the behavior of the system is symmetric. The two oscillators exchange 
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an equal amount of energy amongst themselves per unit time. The interaction Hamiltonian takes both oscillators 
out of their respective equilibria and, eventually, after many iterations of the process, the two attain an equilibrium 
Maxwcll-Boltzmann distribution at an effective temperature Ti^oo) = pi(0) + 72(0)] /2 that is different from their 
initial temperatures. The equilibrium reached here is a thermal equilibrium, in other words, the oscillators have 
reached a state where their temperatures have ceased to change, and a single temperature can be attributed to the 
entire system. This was in fact noted by several individuals including Montroll and Shuler 
and Falkoff 



29], Mathews, Shapiro, 



30| | , as well as Rau [X 61 ] , but the availability of numerical techniques for carrying out the iteration was 
not present at the time. Nonetheless, the predictions made were correct and the results are not surprising. They are 
simply a consequence of the zeroth law of thermodynamics: when two systems are put in contact with each other, 
there will be a net exchange of energy between them unless or until they are in thermal equilibrium. This phenomenon 
has also been studied by Andersen and Shuler for the specific case of the relaxation of a hard-sphere Rayleigh and 
Lorentz gas [31 1. 



In all of the above examples, the frequencies of the oscillators were equal. A somewhat more unexpected result 
(from a thcrmodynamical point of view), appears in the case where u)\ ^ u)%. Here, the exchange of energy in the 
system is not symmetric, such behavior being caused by the fact that the frequencies are not equal. The equilibrium 
reached is not a thermal one and the final combined distribution no longer satisfies Maxwell-Boltzmann statistics. 
The oscillators interact through the interaction Hamiltonian which takes them out of their initial equilibrium, they 
exchange quantities (including temperature), and eventually reach a state where each oscillator has attained a new 
equilibrium at it's own temperature, such that w\/\kT\{po)] = (^/[fcT^oo)]. Besides imposing the condition that the 
total number of particles is a constant, the interaction Hamiltonian also dictates that the free energy of the system 
(the energy not including the interaction energy) is conserved only when the frequencies of the oscillators are equal. 
Here, the selection rules imposed by the interaction Hamiltonian override the statistical mechanical effects. 



V. CONCLUSIONS AND BEYOND 



We have explicitly calculated the time evolution of an initially uncoupled system of two harmonic oscillators, under 
the repeated application of an interaction Hamiltonian for successive time intervals r. Our results were calculated 
exactly using iterative methods and without using perturbation theory. We started with two initially uncoupled 
quantum harmonic oscillators, each at equilibrium and with temperatures T\ and Ti respectively, and showed that 
after enough periodically repeated applications of the interaction Hamiltonian, the oscillators came to equilibrium 
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with each other - if their frequencies were equal - at an effective temperature which was different from their initial 
temperatures. When the frequencies were unequal, each oscillator attained a new equilibrium but there was no 
effective temperature that could be attributed to both of them. On the one hand, when the frequencies of the 
oscillators were the same, the system's initial Maxwell-Boltzmann distributions evolved into a thermal equilibrium 
at a new Maxwell-Boltzmann distribution through a series of transient Boltzmann distributions and on the other 
hand, when the oscillator frequencies were unequal, the equilibrium reached was not thermal and the final combined 
distribution at equilibrium no longer satisfied Maxwell-Boltzmann statistics. To explain this, we showed that, at 
equilibrium, wi#i(oo) = ^26*2(00). Further, we have concluded that the selection rules imposed by the interaction 
Hamiltonian override statistical mechanical effects. 

An important ingredient of our model is what we have called the "refresh" procedure. After applying a constant 
interaction Hamiltonian to our initially uncoupled set of harmonic oscillators, we perform an act of "refreshing". This 
corresponds to the assumption that there are no fixed phase relationships between the two harmonic oscillators, while 
at the same time, the phase relationships between the number states of each individual oscillator are definite. In 
addition, we assume that we have no knowledge about the occupation numbers of each of the modes in the second 
system and so we average over all these occupation numbers. The fact that the quantity N = n\ + n 2 is conserved 
by definition of the total Hamiltonian in our model, and the fact that the initial density matrix is diagonal in the 
number basis {\nin 2 )} of the harmonic oscillators, forces the evolved reduced density matrices /9i(nr) and p 2 (nr) 
with n = 0, 1, 2, 00, to be diagonal in their respective number bases {|^i} and {|«2}- In turn, this allows us to 
write the new combined density matrix after the "refreshing" as a product state of the new reduced states of the two 
oscillators. We have shown how the "refresh" procedure can be obtained through a single, non-selective measurement 
on one of the two oscillators. A single spectrum measurement in the number basis of oscillator 1, results in the 
diagonalization of the output state in both oscillator number bases {|ni)} and {|n2)}, and this in effect permits us to 
assume statistical independence between the two post-measurement states by allowing us to write the global output 
state as a tensor product of pi and p 2 . Even though the transformation of the "refresh" procedure formally appears 
to be nonlinear in pi 2 , this is simply a consequence of how we choose to express it. In effect, the "refreshing" model 
in no different than performing a spin measurement on an entangled electron pair. Naturally, the results obtained in 
this paper are dependent on the specific model presented, and more specifically, the interaction Hamiltonian and the 
"refresh" procedure. The question we may now ask is: How will our results change if any or both of these factors are 
changed? 
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